MCNP Deck Generator

Generates multiple decks to be ran in MCNP for the purposes of determining the sensitivity of the CCD+lanex used in the electron spectrometer. All decks include the lanex layers and use a monoenergetic source of electrons. Each deck varies in the electron energy used as well as the angle of incidence


In [2]:
#Imports
from math import *
import cmath
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import sys
import os

#Import custom modules
from physics import *

In [110]:
#Defined functions
B0 = 1710.0 #Magnetic field strength in Gauss
B0 = B0/10**4
yM = 0.5*.0254

def Radius(KE):
    '''Radius of electron orbit in m given KE in keV'''
    return me*c/(q*B0)*sqrt(1-(1/(KE*1000*q/(me*c**2)+1))**2)*(KE*1000*q/(me*c**2)+1)

def AngleIncidence(KE):
    R = Radius(KE)
    return asin((R-yM)/R)

In [116]:
#Create arrays over 2D parameter space (KE and AOI)
KEstepsize = 0.5
KE = np.arange(0.1,3.0+KEstepsize,KEstepsize)
ThetaStepsize = 30.0
ThetaStepsize *= 2*pi/360
Theta = np.arange(0,pi/2+ThetaStepsize,ThetaStepsize)
CosTheta = np.round(np.cos(Theta),6)

Nsegments = 100
SegmentPosition = np.linspace(0.0255,0.0365,Nsegments+2)
SegmentingSurfaces = ''
SurfaceSegmentNumber = range(100,100+Nsegments)
SurfaceSegmentNumberString = ''
line = ''
for x in range(Nsegments):
    if len(line)<=45:
        line += '-' + str(SurfaceSegmentNumber[x]) + ' '
        SurfaceSegmentNumberString += '-' + str(SurfaceSegmentNumber[x]) + ' '
    else:
        line = ''
        SurfaceSegmentNumberString += '\n\t-' + str(SurfaceSegmentNumber[x]) + ' '
    
    SegmentingSurfaces += '   {Segment}       pz {Position}\n'.format(Segment=100+x,Position=SegmentPosition[x+1])

for i in KE:
    for j in CosTheta:
        directory = os.curdir + '/MCNP_Decks/MCNP_{KE}MeV_{Theta}Degrees'\
        .format(KE=str(round(i,3)),Theta=str(round(np.arccos(j)*360/(2*pi),2)))
        file = open(directory, "w")
        file.write('''c  ************************
c         Cell Cards
c  ************************
c  number  material density then logical definition with surfaces
    1   1   -2.70  1 -2 6 -7 8 -9 $Al layer
    2   2   -1.44  2 -3 6 -7 8 -9 $Cellophane layer
    3   3   -1.38  3 -4 6 -7 8 -9 $PET layer
    4   4   -4.48  4 -5 6 -7 8 -9 $Phosphor layer
    5     5  -0.000020615 -999 #1 #2 #3 #4
    6     0         999 

c  ************************
c       Surface Cards
c  ************************
c number    type  position 
    1        pz 0 $Lanex surface
    2        pz 0.0025 $End of Al layer
    3        pz 0.008 $End of cellophane layer
    4        pz 0.0255 $End of PET layer
    5        pz 0.0365 $End of phosphor layer
    6        px -0.1 $Lanex left edge
    7        px 0.1 $Lanex right edge
    8        py -1.5 $Lanex bottom edge
    9        py 1.5 $Lanex top edge\n''')
        file.write(SegmentingSurfaces)
        file.write('''   999       so 100

c  ************************
c        Problem Type
c  ************************
mode p e
c  **********************************
c        Material Specification
c  **********************************
m1  13000.02p -1  $aluminum  -2.70
m2  1001.03e -0.062162  $cellophane  -1.44
     6012.03e -0.444462
     8016.03e -0.493376
m3  1001.03e  -0.143711  $polyethylene -.9
      6000.03e  -0.856289 
m4  64000.03e 2  $Phosphor  -4.48
     8016.03e 2
     16000.03e 1
m5  1001.03e  2  $water vapor  -0.000020615
     8016.03e  1 
c
imp:p             1 4r                      0  $ 1, 3
c
c  *********************************
c       Physics Simplification
c  *********************************
phys:p  100 0 0 .9 0 $ pair production on change -1 to # for bias    
phys:e  100 0 0 1 0 20 0 1 1 0   $ bremstrahlung biased  
c  *************************
c           Tally 
c  *************************
F11:e 5   $particle counting
F16:e 4   $Energy deposition averaged over a cell\n''')
        file.write('FS16 {Segments}T\n'.format(Segments=SurfaceSegmentNumberString))
        file.write('''c  **************************
c      Energy Tally Card 
c  **************************
c e0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2
c     1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4
c     2.5 2.6 2.7 2.8 2.9 3.0
c  ************************
c     Source Definition
c  ************************\n''')
        file.write('sdef DIR={mu} VEC=0 0 1 POS=0 0 -0.1 CEL=5 ERG={KE} PAR=3\n'.format(KE=str(i),mu=str(j)))
        file.write('''c sp1 -5 1
c
c Number of histories to transport
nps 20000\n''')

        file.close()

In [139]:
#Create arrays over 1D parameter space (KE and AOI associated with particle tracking simulations)
KEstepsize = 0.01
KEsim = np.arange(1.00,3.00+KEstepsize,KEstepsize)
CosThetaSim = []
for x in KEsim:
    CosThetaSim.append(np.cos(AngleIncidence(x*1000)))

Nsegments = 100
SegmentPosition = np.linspace(0.0255,0.0365,Nsegments+2)
SegmentingSurfaces = ''
SurfaceSegmentNumber = range(100,100+Nsegments)
SurfaceSegmentNumberString = ''
line = ''
for x in range(Nsegments):
    if len(line)<=45:
        line += '-' + str(SurfaceSegmentNumber[x]) + ' '
        SurfaceSegmentNumberString += '-' + str(SurfaceSegmentNumber[x]) + ' '
    else:
        line = ''
        SurfaceSegmentNumberString += '\n\t-' + str(SurfaceSegmentNumber[x]) + ' '
    
    SegmentingSurfaces += '   {Segment}       pz {Position}\n'.format(Segment=100+x,Position=SegmentPosition[x+1])

for i in range(len(KEsim)):
    directory = os.curdir + '/MCNP_Decks/MCNP_{KE}MeV_{Theta}Degrees'\
    .format(KE=str(round(KEsim[i],3)),Theta=str(round(np.arccos(CosThetaSim[i])*360/(2*pi),1)))
    file = open(directory, "w")
    file.write('''c  ************************
c         Cell Cards
c  ************************
c  number  material density then logical definition with surfaces
    1   1   -2.70  1 -2 6 -7 8 -9 $Al layer
    2   2   -1.44  2 -3 6 -7 8 -9 $Cellophane layer
    3   3   -1.38  3 -4 6 -7 8 -9 $PET layer
    4   4   -4.48  4 -5 6 -7 8 -9 $Phosphor layer
    5     5  -0.000020615 -999 #1 #2 #3 #4
    6     0         999 

c  ************************
c       Surface Cards
c  ************************
c number    type  position 
    1        pz 0 $Lanex surface
    2        pz 0.0025 $End of Al layer
    3        pz 0.008 $End of cellophane layer
    4        pz 0.0255 $End of PET layer
    5        pz 0.0365 $End of phosphor layer
    6        px -0.1 $Lanex left edge
    7        px 0.1 $Lanex right edge
    8        py -1.5 $Lanex bottom edge
    9        py 1.5 $Lanex top edge\n''')
    file.write(SegmentingSurfaces)
    file.write('''   999       so 100

c  ************************
c        Problem Type
c  ************************
mode p e
c  **********************************
c        Material Specification
c  **********************************
m1  13000.02p -1  $aluminum  -2.70
m2  1001.03e -0.062162  $cellophane  -1.44
     6012.03e -0.444462
     8016.03e -0.493376
m3  1001.03e  -0.143711  $polyethylene -.9
      6000.03e  -0.856289 
m4  64000.03e 2  $Phosphor  -4.48
     8016.03e 2
     16000.03e 1
m5  1001.03e  2  $water vapor  -0.000020615
     8016.03e  1 
c
imp:p             1 4r                      0  $ 1, 3
c
c  *********************************
c       Physics Simplification
c  *********************************
phys:p  100 0 0 .9 0 $ pair production on change -1 to # for bias    
phys:e  100 0 0 1 0 20 0 1 1 0   $ bremstrahlung biased  
c  *************************
c           Tally 
c  *************************
F11:e 5   $particle counting
F16:e 4   $Energy deposition averaged over a cell\n''')
    file.write('FS16 {Segments}T\n'.format(Segments=SurfaceSegmentNumberString))
    file.write('''c  **************************
c      Energy Tally Card 
c  **************************
c e0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2
c     1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4
c     2.5 2.6 2.7 2.8 2.9 3.0
c  ************************
c     Source Definition
c  ************************\n''')
    file.write('sdef DIR={mu} VEC=0 0 1 POS=0 0 -0.1 CEL=5 ERG={KE} PAR=3\n'\
               .format(KE=str(KEsim[i]),mu=str(CosThetaSim[i])))
    file.write('''c sp1 -5 1
c
c Number of histories to transport
nps 20000\n''')

    file.close()